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dimensional quantum dots is presented. Advantages and limitations are studied through comparison 
with other high accuracy approaches for two to eight confined electrons. The possibility to effectively 
use a very large basis set is found to be an important advantage compared to full configuration in- 
teraction implementations. For the two to eight electron ground states, with a confinement strength 
close to what is used in experiments, the error in the energy introduced by truncating triple exci- 
tations and beyond is shown to be on the same level or less than the differences in energy given by 
two different Quantum Monte Carlo methods. Convergence of the iterative solution of the coupled 
cluster equations is, for some cases, found for surprisingly weak confinement strengths even when 
starting from a non-interacting basis. The limit where the missing triple and higher excitations 
become relevant is investigated through comparison with full Configuration Interaction results. 
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I. INTRODUCTION 

Ever since Tarucha et al^ experimentally demon- 
strated atom-like properties in few electron quantum 
dots, in particular the existence of a shell structure, these 
systems have attracted a lot of theoretical interest as new 
targets for many-body methods. In contrast to the sit- 
uation for the naturally occurring many-body systems, 
the strength of the overall confinement relative to that 
of the inter-particle interaction can here be freely varied 
over a large range, at least in theory, and thus completely 
new regimes can be explored. When the aim is to study 
the performance of a specific many-body method, it is 
justified to use a simple model for the confining poten- 
tial and most theoretical studies restrict themselves to 
the two dimensional harmonic oscillator potential. The 
interaction between the dot electrons and the surround- 
ing semi-conductor material is further usually modeled 
through the use of a material specific effective electron 
mass and relative dielectric constant. This model implies 
thus a two-dimensional truly atom-like device, on which 
calculational methods developed for atoms can be ap- 
plied after minor adjustments. Early calculations proved 
the model to be adequate. Combined with a reason- 
able account for the electron-electron interaction through 
methods such as Density Functional Theory (DFT)^"— or 
Hartree-Focki^— , the two-dimensional harmonic oscilla- 
tor confining potential did indeed give a good qualitative 
agreement with experiments. In particular the closed 
shells forming with two, six and twelve trapped electrons 
could be explained, as shown in many studies, see e.g. 
the review by Reimann and Mannincn^^. 

Neither the true form of the dot confinement, nor the 
extent to which it deviates from being purely two di- 
mensional, is easily extracted from experiments. If the 
electron-electron interaction could be sufficiently well ac- 
counted for, it might be possible to gain further under- 
standing of such properties through comparison between 



experiment and theory. Some studies in this direction 
have been performed, e.g. by Matagne et al}"^ who made 
quantitative statements about the non-harmonicity of the 
confining potential from the comparison between DFT 
calculations and experiments. DFT has proven to be 
able to account for a large part of the electron corre- 
lation, but still it is not really the best choice for such 
investigations since it is hard to make an a priori esti- 
mate of the obtained error. There are several ways to 
account for correlation more systematically. A number 
of studies on quantum dots and related structures have 
been carried out with Configuration Interaction (CI), see 
e.g. RefsJ^"— . Full CI is in principle exact and ap- 
plicable for all relative interaction strengths. The term 
"full CI" refers to a calculation where all Slater deter- 
minants, obtained by exciting all possible electrons to 
all possible orbitals that are unoccupied in the studied 
electronic configuration, are included. It is obvious that 
the size of the full CI problem grows extremely fast both 
with the number of particles and with the size of the 
basis set (used to represent the unoccupied orbitals). It 
is well known that truncated CI, i.e. keeping e.g. only 
single and double excitations from the leading configu- 
ration, lacks size-extensivity. In short this implies that 
the method does not scale properly with the size of the 
studied system, see for example the review by Bartlett^. 
Truncation of the number of excitations is thus not a real 
option, and instead all but the smallest systems have to 
be calculated with very small basis sets. Recently a thor- 
ough investigation of the performance of full CI applied 
to quantum dots with a basis set consisting of harmonic 
oscillator eigenfunctions was made by Kvaali^ and the 
main conclusion was that the convergence with respect 
to the size of the basis set was slow and that additional 
features such as effective two-body interactions have to 
be added for meaningful comparisons with experiments. 
An alternative approach for high accuracy calculations 
are Quantum Monte Carlo (QMC) methods, which sue- 
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cessfully have been applied to quantum dotg--"--. Here 
the computational cost grows modestly with the number 
of electrons, and the method provides a very efficient way 
to calculate the ground state for a specific symmetry. The 
nodal structure of the trial wave function can be used to 
impose restrictions on the solutions so that also excited 
states can be obtained to some extent, see e.g. the discus- 
sion in the review by Foulkes et al.^'^. Still, calculations 
on general excited states are not straightforward and ad- 
ditional methods are needed for realistic calculations on 
important parts of the quantum dot spectrum. 

Most of the methods used on atoms and molecules have 
also been applied to quantum dots in several implementa- 
tions. The least studied method is however many-body 
perturbation theory which has been shown to be very 
powerful for the calculation of atomic properties. Calcu- 
lations up to second order in the perturbation expansion 
have been made by just a few authors^^i^^, and equally 
few coupled-cluster calculations have been presented^^^. 
For small to medium sized molecules, as well as atoms, 
the coupled-cluster(CC) method is known to success- 
fully combine feasibility with accuracy. The coupled- 
cluster theory was introduced in 1960 by Coester and 
Kiimmel-- in nuclear physics, and since then contribu- 
tions have been given by many authors. A rather recent 
review regarding its performance in quantum chemistry 
has been made by Bartlett and Musiai^. We present 
here a thorough investigation of how the coupled-cluster 
method with single and double excitations (CCSD) per- 
forms, in comparison with full CI and Quantum Monte 
Carlo-studies, on two-dimensional quantum dots. In sec- 
tion|lI]we summarize CCSD and briefly discuss its advan- 
tages. In section Hill our implementation for calculations 
on circular quantum dots is outlined. In Section IIVI we 
present results for dots with two to eight electrons and 
compare them with those obtained with other methods. 
The first question is whether the restriction to single and 
double excitations is adequate. It is known to be a good 
approximation for atoms, but we expect it to eventually 
fail for sufficiently weak confinement strengths. Here, we 
try to establish when this happens. The next point is 
the feasibility and we show that results converged with 
respect to e.g. basis size can be obtained for much larger 
systems than is possible for CI calculations. We show for 
hid ^ 3 meV and the = 2 to 8 groundstates that the 
error relative to Diffusion Monte Carlo results is on the 
same level or less than the difference between the energies 
given by the Variational and Diffusion Quantum Monte 
Carlo methods. 



II. THEORY 

The formalism used in the present study can be 
found in more detail in the textbook by Lindgren and 
Morrison^. Here we just discuss the aspects important 
for the understanding of the present results. In order to 



solve the Schrodinger equation 

iJf = E^, (1) 

for an N -ferniion system, the Hamiltonian is partitioned 
as 

H = Ho + V, (2) 

where the eigenstates of Hq are known and V is the re- 
mainder, i.e. it is the perturbation with respect to the 
already solved Hamiltonian Hq. In the present study Hq 
is usually the Hamiltonian for the non-interacting sys- 
tem, and thus V is the whole electron-electron interac- 
tion, but also other choices have been examined as will 
be discussed below. We will further assume that Hq is a 
sum of N single-particle Hamiltonians with known eigen- 
states: 

N 

Hq = ^ hi , 

i 

hi\i) = £i I ?). (3) 

With the orbitals, | i), we can form a model space suitable 
for the state (s) we are interested in. In the following we 
will use a one-dimensional model space, P, spanned by 
one Slater determinant 

a = {abcd...N}, (4) 

P = I a){a I, (5) 

where a,b., c,d . . . , N denote the occupied orbitals, and 
the curly brackets denote antisymmetrization. A multi- 
dimensional extended model space is also a possibility^, 
which though will be left for future investigations. The 
model function is the projection of the exact solution, 
Eq. ([ij, onto the model space 

*o = (6) 

We further assume that the model function, but not the 
full wave function, is normalized, a condition usually re- 
ferred to as intermediate normalization. It is possible 
to define a wave operator, that transforms the model 
function into the exact state, i.e. 4' — 57^o- To separate 
the part of that projects onto the model space and 
that which brings the solution out of the model space we 
write 

^ = l + X, (7) 

where x is sometimes referred to as the correlation op- 
erator. The wave operator can be obtained from the 
generalized^ Block equatiorJ^^, which we will use in the 
form^iil; 

[n,HQ]p = {Qvnp-xPvnp), (8) 

where Q is the orthogonal space to P, such that P + Q = 
1. Eq. ([5]) is equivalent to the Schrodinger equation, but 
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in this form it allows for an iterative solution procedure. 
Setting X on the right hand side initially to zero, one 
can obtain a first approximation of x which then can be 
inserted on the right-hand side to get a better approx- 
imation and so on until convergence is reached. When 
an order by order expansion is carried through, it can be 
shown that in each order the so called unlinked diagrams 
from the first term in Eq.|S]are canceled by contributions 
from the second term. Unlinked contributions are those 
that include parts that are surrounded by P-operators, as 
thexPFriP-term in Eq. [5] This is the so called linked di- 
agram theorem^ see for example Refs<SSi^o and references 
therein. It is thus possible to keep only linked contribu- 
tions in the expansion. While the exclusion principle is 
obeyed for the sum of linked and unlinked contributions, 
the cancellation of the unlinked contributions is achieved 
only if the exclusion principle is lifted, resulting in ex- 
clusion principle violating diagrams also in the retained 
linked contributions. It is the linked expansion which will 
be the basis for the coupled cluster expansion below. 

When Q. has been obtained it can be used to construct 
the effective Hamiltonian 



H^g ^ phqP + pvnp, 



(9) 



which gives the exact energies when acting on the model 
spaced. The total energy can then be written as 



E - {a\H,B\a) = {a\Ho + V + Vx\a), 



(10) 



where P | a) =| a) has been used. 

Q is the complementary space to P and can formally 
be built up from all Slater determinants, /3, that differ 
from a 



(11) 



The space spanned by Q is in principle infinite, but in 
practice we use a finite basis set to represent the eigen- 
states to h, Eq. This makes also the Q-space finite, 
but still it grows rapidly both with the size of the basis 
set and with the number of particles. We focus now on 
the representation of the Q-space for several particles. 
For this purpose we can classify the Slater determinants 
belonging to Q with respect to by how many single par- 
ticles states they differ from P. For example 



al = {rbcd ...N} 



(12) 



differs from Eq. ([4]) only in that a has been replaced by 
r, and it is labeled a single excitation, while 



a2^{rscd...N}., 



(13) 



differs from Eq. ^ in that a and b have been replaced by 
r and s, and it is labeled a double excitation. A complete 
calculation on a two-particle system requires single and 
double excitations, while such a calculation on a three 
particle system also requires triple excitations, and so 



on. For a general many-particle system it is necessary to 
truncate this series at some point due to both complexity 
and computational load. For this truncation there exists 
several choices, x is the part of the wave operator that 
lies in the Q-spa,ce. It can for example be divided up as 



X = Xl + X2 + X3 + 



(14) 



where the subscripts denote the number of excitations. 
If we truncate this sum after e.g. xa, we will reproduce 
CI with single, double, and triple excitations^. With 
the coupled-cluster approach the truncation is made in 
an alternative way. First we start from the linked form 
of the Bloch equation 



[n, Hq] p = {Qvnp - xPvnp) 



linked ' 



(15) 



where only linked contributions are retained in the itera- 
tive procedure. As a second step we define a cluster op- 
erator S — Si-\-S2 + S3-\-..., where each term represents 
the connected part of the wave operator for n excitations, 



connected 



The term connected denotes that the 
wave operator cannot be divided up in parts where the 
particles interact independently in smaller clusters, e.g. 
two-by-two. The S operator can be shown— to satisfy a 
Bloch type equation 



[s, Ho] p = (Qvnp - xPvnp) 



connected " 



(16) 



The wave operator, f2, can now be written as a sum of 
products of ^n-operators. All such terms are generated 
through the exponential ansatz: 

n = {expiS}} = l+Si+S2+^{Siy+{SiS2}+^{Sf} 



1 



1 



{Si} + ^{s!s2} + ^{sf} 



4! 



(17) 



The curly brackets denote here that it is the normal or- 
dered products of the operators that should be used, 
which implies antisymmetrization. We can now identify 
all single, double, triple etc excitations accordingly 



^2 



Si 



S2 + ^{S'i} 



S4 -\- {S1S3} 



l{si} + l{s!s2} + ^{st} 



(18) 



From Eq. [TH] it is clear that when truncating after the 
S2 cluster, we still include the parts of and il^ that 
can be written as combinations of Si and ^2 operators, 
i.e. the terms in boxes above. This is the Couple Cluster 
Singles and Doubles method. How these products of S- 
operators enter in the expansion will become more clear 
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when Eqs. (I29ll30p are discussed below. See also Ref.^ for 
more details. There are two clear advantages of this trun- 
cation scheme. First, the probably most important triple 
and quadruple excitations are now included in a scheme 
that is much less computationally demanding than calcu- 
lating full triples and quadruples. Second, and this is in 
contrast to the scheme indicated in Eq. (|14p , the inclusion 
of the disconnected products makes the coupled-cluster 
method size extensive also in its truncated version. 

In the following we will investigate the performance of 
the coupled cluster method when including all Si and 
52 terms in Eg. dTSl) (the expressions in boxes), i.e. the 
CCSD method. Although the practical implementation 
is different, the present study includes the same effects 
as the implementation for atoms by Salomonson and 
Oster'*'*, where more details about the method can be 
found. 



III. IMPLEMENTATION 



order is 10 and combined with the knot sequence this 
yields 29 radial basis functions, UnmimA''^)^ fo'" each com- 
bination (rni,ms). The lower energy basis functions are 
physical states, while the higher ones are mainly deter- 
mined by the box. The unphysical higher energy states 
are, however, still essential for the completeness of the 
basis set. 

Eq. (|2m22p implies that the one-particle Schrodinger 
equation, (jl9p . can be written as a matrix equation 



he = eBc 



(23) 



where /i^ = (B,e™n^l-S. 



and B,j = 



Eq. I|23p is a generalized eigenvalue problem that can 
be solved with standard numerical routines. The inte- 
grals in (1231) are calculated with Gaussian quadrature. 
B-splines are piecewise polynomials, and since also the 
potential is in polynomial form in Eq. (|19p . essentially 
no numerical error is produced in the integration. 



A. Single-particle treatment 

For a single particle confined in a circularly symmetric 
potential the Hamiltonian reads 



2m* 



(19) 



where the effective electron mass is denoted with m*. 
With a pure harmonic confinement we have 



B. The Coulomb interaction 

The perturbation V in Eq. ([2]) will include the electron- 
electron interaction not accounted for in Hq. It can be 
the full Coulomb interaction or the difference between 
that and some mean-field approximation. In either case, 
we need a suitable way of dealing with the Coulomb inter- 
action in two dimensions. As suggested by Cohl et al^, 
the inverse radial distance can be expanded in cylindrical 
coordinates {R, (/>, z) as 



Uc{r) = —m*uj^r^. 



(20) 



This is the confining potential used in all numerical re- 
sults in the present study but any circularly symmetric 
confinement can be used in the developed computer code. 

The single particle wave functions separate in polar 
coordinates as 



( r, ( 



(r)e 



imi <p I, 



(21) 



We expand the radial part of the wave functions in so 
called B-splines labeled Bi with coefficients Ci, i.e. 



..(0 =E^*^■'W■ 



(22) 



1=1 



B-splines are piecewise polynomials of a chosen order fc, 
defined on a so called knot sequence and they form a com- 
plete set for the linear space defined by the knot sequence 
and the polynomial order*^. Here we have used 40 points 
in the knot sequence, distributed by the use of an arcsin- 
function. The last knot point, defining the boundary of 
the box to which we limit our problem, is scaled with the 
potenti al strengt h through the harmonic oscillator length 
unit yjh/ {m*uj). For example with hiu sa 11.857 meV 
(which corresponds 1 a.u.* for GaAs, see Section ITV)) the 
last knot point is located at r ~ 70 nm. The polynomial 



1 



1 



|ri-r2| TT^/WRa 



im((l>i-4i2) 



where 



X 



Rl + Rl + {zi- 
2R1R2 



(24) 



(25) 



Assuming a two-dimensional confinement we set zi = z^ 
in (pSj). The Q,„_ 1 (x)-functions are Legendre functions 
of the second kind and half-integer degree. We evaluate 
them using a modified^ version of software DTORHl.f 
described in^^. 

Using Eqs. (PT|) and (1^ . we can write the electron- 
electron interaction matrix element as 

{ah\^\cd) = 

{ua {ri)ub irj)\ 1^'' \uc{n)udirj)) 



E 



X (m^|m^)(™t|™f). 



(26) 



Note that the angular (0) integration in (|26p yields a 



non-zero result only if m 



ran 



ran 



ran. This 
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determines the degree m of the Legendre-function in the 
radial part of Eq. (f^ . It is also clear from Eq. (pS)) 
that the electron-electron matrix element equals zero if 
orbital a and c or orbitals b and d have different spin 
directions. 



C. Alternative single particle potentials 

If V in Eq. ^ has to account for the full electron- 
electron interaction, as it does when the single parti- 
cle Hamiltonian just includes the external confinement 
potential, Eq. (|20|) . it might be difficult to obtain con- 
vergence of the iterative solutions of Eq. ([8]). At least 
this is expected for weak external confinement and for 
many-electron dots. A remedy is then to start from a 
Hamiltonian Hq that already includes the bulk of the 
electron-electron interaction. Many choices are possible 
here and we have investigated two of these. For a sin- 
gle Slater determinant the Hartree-Fock approximation 
is known to minimize the energy. In this approxima- 
tion, each electron moves in the average potential from 
the other electrons. The one-particle potential, including 
the external confinement, is then 



.HF 



= {Bj\u,{r)\B, 



+ 



a<N 



[B,a\^\B,a) ^ {B,a\^\aB,). (27) 
ri2 ri2 



The last term gives the (non-local) exchange interaction. 
One complication with Eq. ([ST)) is that for a situation 
where not all electron spins are paired, electrons with 
the same quantum numbers n and mg but with different 
spin directions will experience different potentials. As 
a consequence, the total spin = (X^i^O^: does not 
commute with the Hartree-Fock Hamiltonian. In spite of 
the fact that the full Hamiltonian, Eq.©, still commutes 
with S^, this property might lead to complications, see 
Refj2^ for more details, and it might be more practical 
to use a starting point where the exchange interaction is 
approximated with a local potential. We adopt for this 
a traditional Local Density Approximation (LDA) with 
a variable amount of exchange 



,LDA 



{Bj\u,{r)\Bi 



-j.<N 



[B,a\^\B,a} 
ri2 



v{Bj\4a*B 



Mr) 



B, 



(28) 



where p{r) is the radial electron density and 77 is the so 
called Slaters exchange parameter, which often is set to 
one. In Section HVl we present results with r; = 1.0 and 
■q = 1.4. 



D. Many-Body treatment 

Equipped with a finite representation of the Q space it 
is possible to construct the ^n-operators, and thus also 



the wave operator fi. We now use the Coupled Clus- 
ter Singles and Doubles truncation of the possible ex- 
citations, i.e. only the terms in the boxes in Eq. (fT8|) 
are kept. We start from Eq. (|16p and note that for a 
model space built from a single Slater determinant, the 
xPVftP-teim is fully canceled by the unlinked diagrams 
from the QFfiP-term. Only the QVilP-teTm remains 
thus on the right-hand side of Eq. ([TB| . Starting with 
57^^) = 1 and x^^^ — 0, we can write the recursion rela- 
tion for the the S'l-amplitudes as 



«\S,\a) 



1 



(a^\Vi + VSi+VS2 



+ ^V{Sf} + V2{S,S2} + l^V2{Sf}\ay (29) 



and for the S'2-amplitudes 



{a:i\S2\a) 



i+l 



1 



+ — — 



C^Z\V2^V2Si+VS2 



1 



1 



1 



2! 

r^2{5ni«) 



(30) 



, ^V2{Sl} + V{S,S2} + ^V2{Sl} + -Msl} 

+ ^V2{SlS2} + i 

where only connected contributions should be kept on 
the right-hand side. Here F = T^i -I- V2 is the total per- 
turbation, Vi is the part of the perturbation that can be 
written as a one-particle operator and V2 is the part of 
the perturbation that can be written as a two particle 
operator. The index i denotes the iteration number. It 
is related, but not equal, to the order in the perturbation 
expansion. The quoted figures in section IIVI are always 
self-consistent with respect to Eqs. (|29ll30| . Note that, 
e.g. the single excitation cluster, 5*1 ([^ . is built from 
up single, double and triple excitations. As an example 
we note that the included triples are those that can be 
be written as disconnected singles connected by a pertur- 
bation V2 (the last term on the second line of Eq. (|29|)). 
as well as combinations of singles and doubles connected 
by V2 (second term on the last line of Eg.ip^V These 
are so called intermediate triple excitations. In a similar 
way also intermediate triples and quadruples contribute 

to 5*2. 

Finally, it is appropriate to comment on the difference 
between the two-dimensional many-body procedure and 
the more studied three-dimensional case, especially with 
respect to the angular integration. The angular momen- 
tum algebra is considerably simplified in two dimensions 
compared to three dimensions. In two dimensions an or- 
bital is defined by only three quantum numbers. With 
polar coordinates these are the radial quantum number, 
n, the angular quantum number, mi, and the spin di- 
rection, mg. The radial functions Unmfm, (?"), Ea.(pT|). 
depend on two of these quantum numbers, n,mi, while 
an additional dependence on rus only arises in case an 
external magnetic field is applied to the dot. In three 
dimensions the desired total angular momentum has to 
be constructed through a linear combination of the dif- 
ferent magnetic components of the orbitals. In spite of 
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the advanced formalisms (e.g. Racah algebra) developed 
in order to avoid explicit summation over magnetic sub- 
states, the angular integration often gets rather cumber- 
some, at least for general open shell configurations. In 
two dimensions there is only place for one particle in each 
spatial orbital and any state with maximized total spin 
can be treated as a closed shell configuration is handled 
in there dimensions. 



IV. RESULTS 

In our numerical studies we use m* = 0.067me and 
er = 12.4 corresponding to the bulk value in GaAs. 



A. Validation 

Table |I] compares the present results with those ob- 
tained by Full Configuration Interaction, Refsj^'^^'^, 
for 2-6 and 8 electrons when starting from the non- 
interacting basis. The purpose of Table |T] is first to com- 
pare with calculations which include exactly the same 
physical effects. Such a comparison can in principle only 
be done for two electrons due to the truncation of S3- 
clusters and beyond in the CCSD method. A second 
purpose is to compare the accuracy and basis set conver- 
gence for more than two electrons and for confinement 
strengths close to the region of interest. 

The key parameter here is the strength of the electron- 
electron interaction relative the confinement provided by 
the harmonic oscillator potential, which can be quanti- 
fied, as in Table U by the length parameter A, 

A=At^ = A^- (31) 
ojm* 47reo£r'i com* 

where m* is the effective mass, is the relative dielectric 
constant, and Og is the effective Bohr radius. We recall 

here that \/ is the harmonic oscillator length unit. 

Larger A-values correspond to a weaker confinement and 
an increased relative importance of the electron-electron 
interaction. For Table U we (mainly) choose A = 2 corre- 
sponding to ftw « 2.964 meV (for GaAs parameters) since 
the next available Cl-results are either much stronger or 
much weaker than this. 

For two electrons, N = 2, both the CCSD and the FCI 
method take all electron-electron effects into account and 
the accuracy is thus only limited by the size of the ba- 
sis set and the numerical procedure. When comparing 
the N — 2 results produced with identical basis sets in 
Table HI we note that our values differ from those pro- 
duced by Kvaali^ at most in the seventh digit, while 
those by Rontani et al}"^ differ in the fifth digit. The 
computer code developed by KvaaiiS, is bench-marked 
to machine precision with exact results and it is reason- 
able to believe that its numerical accuracy is the highest. 



The leading numerical errors in the present implemen- 
tation are due to the precision in the (5m-i/2"functions 
produced by DTORHl.f described in^, and their inte- 
gration in Eq. ipSl) . For more than two electrons these 
numerical errors are much smaller than the errors intro- 
duced through truncations, e.g. of basis sets or of S3 
clusters and beyond, and are of no significance. For the 
the remaining part of Table U we restrict the display of 
our results to six digits. 

Since coupled-cluster is an iterative and perturbative 
method, convergence is never guaranteed. We note that 
for iV = 2 in Table [H convergence is still obtained for 
potential strengths as low as A = 8 even though the full 
electron-electron interaction is here taken as the pertur- 
bation. We emphasize that this constitutes a truly non- 
perturbative case; the total energy is almost seven times 
as large as the strength of the confining potential. 

Continuing to > 2, we conclude that the largest 
relative error, calculated as the percentage of the total 
energy, arises for = 3. Still the deviation is never 
larger than 1.5 x 10~^ in units of or ~ 0.1 percent of 
the total energy. For 6 and 8 electrons we also increased 
the basis set significantly beyond what so far has been 
feasible with FCI. It is clear, at least for N > 6 and 
the potential strengths studied here, that the error made 
by the truncation of the Coupled Cluster expansion to 
include only the Si and S2 cluster operators, see Eq. p^ . 
is far smaller than the error made by truncating the basis 
set in the CTcalculations. 

Table U shows only results obtained with the whole 
electron-electron interaction as a perturbation. The rea- 
son for this is the wish to be able to compare with CT 
calculations. These typically use a non- interacting basis 
set, which further is severely truncated. With this start- 
ing point, convergence of the coupled-cluster expansion 
for weaker confinements than A = 2, for iV > 2 can be 
problematic. However, convergence can be obtained for a 
much wider range of confinement strengths with a better 
starting point: e.g. Hartree-Fock or Local density, but 
fair comparison can then only be made with converged 
results. 

A strict limitation with the CCSD approach is the ne- 
glect of true triples, S3 clusters, and beyond. For suffi- 
ciently weak confinements this approximation will dom- 
inate the error. Tabic [Til shows a comparison between 
the CCSD and FCI methods for three and six electrons 
and for a large range of confinement strengths. The pur- 
pose is here to establish how important the limitation to 
Si and 5*2 clusters is. Confinement strengths as weak as 
possible, but still leading to a converged coupled-cluster 
expansion with a non-interacting basis have been used. 
We note that for A < 1 the CCSD-method yields re- 
sults accurate enough for most practical purposes. For 
A = 2 the error due to the neglected effects in the CCSD- 
mcthod is still so small that (keeping table |T] in mind) the 
possibility to use larger basis sets than in a CI calcula- 
tion well compensates for the lack of triples and beyond. 
For = 3 and A = 4 the CCSD did not converge for 
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TABLE I: Comparison between the Coupled Cluster Singles and Doubles method starting from the pure one-electron basis 
(this work) and Full Configuration Interaction according to the software developed by KvaaliS and according to Rontani et 
al}'^ . For truncation of the basis set the so called shell truncation parameter R — 2n + \mi \ is used. Energies are given in units 



of hoj and the number of confined electrons is 


2-6 and 8. A is 


defined in Ea.(PT|l. 














CCSD 




FCI 


125* Ml) 


Confinement Strength 


fiaj(meV) 


Basis set 


This work 


Kvaal 


Rontani et al. 


N=2 |00> 


A = 1.0 


11.85720 


R=5 


3.013625 


3.013626" 










R=6 


3.011019 


3.011020" 










R=7 


3.009234 


3.009236" 






A = 2.0 


2.964301 


R=5 


3.733597 


3.733598" 


3.7338 








R=6 


3.731057 


3.731057" 


3.7312 








R=7 


3.729323 


3.729324" 


3.7295 


|21) 


A = 2.0 




R=5 


4.143592 


4.143592" 


4.1437 








R=6 


4.142946 


4.142946" 


4.1431 








R=7 


4.142581 


4.142581" 


4.1427 


|00) 


A =6.0 


0.3293668 


R=5 


5.784651 




5.7850 


100) 


A =8.0 


0.1852688 


R=5 


6.618102 




6.6185 








R=6 


6.618091 




6.6185 








R=7 


6.618089 




6.6185 


N=3 111) 


A =1.0 


11.85720 


R=6 


6.37600 


6.374293'' 










R=7 


6.37293 


6.371059'' 










R=8 


6.37069 


6.368708'' 










R=10 


6.36773 


6.365615'' 






A =2.0 


2.964301 


R=5 


8.18306 


8.175035'' 


8.1755 








R=6 


8.17896 


8.169913'' 










R=7 


8.17635 


8.166708'' 


8.1671 


N=4 |20) 


A =2.0 


2.964301 


R=7 


13.635 




13.626 


N=5 111) 


A =2.0 


2.964301 


R=5 


20.3697 




20.36 








R=6 


20.3554 




20.34 








R=7 


20.3467 




20.33 


N=6 |00) 


A =2.0 


2.964301 


R=5 


28.0161 


28.0330' 


28.03 








R=6 


27.9912 












R=7 


27.9751 




27.98 








R=10 


27.9529 












R=15 


27.9390 






N=8 |20) 


A =2.0 


2.964301 


R=5 


47.13801 




47.14 








R=10 


46.70369 












R=15 


46.67960 







"Simcn Kvaal 

'Simen Kvaal****, obtained with the software in Ref™ 
■^Patrick Mcrlof*'^ using the software in Ref™. 



the 125*^/^1) = |11) ground state. For this weak con- 
finement the first excited state is still reproduced well. 
Generally we see that this first excited state, which is 
not as localized as the ground state, is reproduced better 
than the ground state throughout the list of confinements 
strengths in Table |TT] . Intuitively this makes sense; true 
triple excitations (and all many-body effects) should be 
relatively more important for more localized states. 



B. Convergence of the basis and the use of 
different starting points 

Fig. [1] depicts the basis set convergence for the N — 2 
ground state using different starting points. The poten- 
tial strength hui = 3.32 meV, corresponding to A « 1.89, 
is chosen to enable comparison with the Quantum Monte 
Carlo results from Ref»22, where it is argued that this 
value is close to the actual values in the experiment by 
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TABLE II: The importance of 5*3 clusters and beyond. The present Coupled Cluster Singles and Doubles results are compared 
to Full Configuration Interactio»i^>^^^ results obtained with the same basis sets. For three electrons the basis is for both 
methods truncated a.t R — 2n + |m;| = 7, and with six electrons it is truncated at 7? = 5. These basis sets are not saturated, 
but the comparison unveils the level at which contributions beyond CCSD contribute, as function of the confinement strength. 
The values in parenthesis are the differences to the corresponding Full Cl-value. Energies are given in units of huj. 



\2S Ml) 


A 


hjj{meV) 


CCSD (present) 


Full CI 


N=3 


111) 

1 / 


0.5 


47.42881 


5.28660 (+0.00019) 


5.28640 






0.75 


21.07947 


5.85048 (+0.00077) 


5.84971 






1.0 


11.85720 


6.37292 (+0.00187) 


6.37106 






1.5 


5.269868 


7.32169 (+0.00539) 


7.31630 






2.0 


2.964300 


8.17635 (+0.00965) 


8.16670 






4.0 


0.7410752 


diverges 


11.0425 




130) 


0.5 


47.42881 


5.90815 (+5x10 ) 


5.90814 






0.75 


21.07947 


6.34019 (+0.00002) 


6.34017 






1.0 


11.85720 


6.75908 (+0.00006) 


6.75903 






1.5 


5.269868 


7.56147 (+0.00020) 


7.56128 






4.0 


0.7410752 


11.0514 ( -0.00120) 


11.0526 


N=6 


100) 


0.1 


1185.720 


11.1979 (+8x10-^) 


11.1979 






0.5 


47.42881 


15.5624 (+0.00062) 


15.5618 






1.0 


11.85720 


20.2609 (+0.00371) 


20.2572 






2.0 


2.964301 


28.0161 ( -0.01687) 


28.0330 



Tarucha et ali^. Fig. [T] shows that the basis sets begin 
to saturate around niax(n) = 10 and for the angular 
quantum number they saturate to the same extent at 
approximately |m;| = 4. The shell truncation param- 
eter R = 2n + |m;|, used in most of the previous CI- 
calculations, is thus here {N = 2 and hui = 3.32 meV) 
clearly overemphasizing the need for high | to/ |- values in 
the basis set, on the expense of n-values. To put it in 
other words; it is apparent from the figure that a ba- 
sis cut of (n, \mi\) < (10,4) is a better choice than e.g. 
R = 2n + I TO; I < 10. For weaker confinements and more 
confined particles the behavior is probably different; one 
expects then the high |to/| basis functions to be relatively 
more important. 

Moreover, Fig. [1] shows that for small n the different 
starting points give quite a spread in the energy but for 
n> 10 the results are virtually independent of the start- 
ing point. Again this demonstrates that for the com- 
parison with Cl-results from KefsM^^^^, which all are 
limited to very few n:s, we are limited to the use of the 
non-interacting starting point. This is unfortunate since 
the CCSD method can handle much weaker confinement 
strengths with a Hartree-Fock or a Local Density Start- 
ing point. 

It is easy to be mislead by Fig. [1] and draw the con- 
clusion that there is no need for other starting points 
than the pure harmonic oscillator basis and that Hartree- 
Fock is the worst of the tested starting points. However, 
the number of needed iterations in Eq. ([29l - [30|) (^ orders 
in the perturbation expansion) to obtain self-consistency 
can differ significantly between the starting points, with 



Hartree-Fock often being the fastest of them all. For ex- 
ample, with a basis size of (n, |to/|) < (28, 1) and includ- 
ing up to second order corrections to the energy, the pure 
harmonic oscillator starting point yields 0.9223 a.u.*, 
while Hartree-Fock gives 1.0449 a.u.* and the fully con- 
verged result is 1.0250 a.u.* (for both starting points). 
That is, with the Hartree-Fock basis we start much closer 
to the fully correlated situation. 

1. Extrapolation and convergence of the basis 

Even though much larger basis sets can be used with 
CCSD than with CI, the calculation time still grows with 
the number of particles and with the size of the basis set. 
Extrapolation of the results to infinitely large basis sets is 
then often an efficient strategy. Many elaborate strate- 
gies can be envisaged here, e.g. adjusting the size of 
the basis set during the iterations. One can for example 
obtain convergence in a limited basis and then systemat- 
ically increase it, or one can filter the contributions and 
only keep those that are estimated to contribute over a 
certain level. Here we restrict ourselves to a brief discus- 
sion of the potential gain of extrapolation. 

A first example is shown in Table IIIIl where the full 
radial basis is used to investigate the rm expansion. The 
extrapolated values are obtained through a linear regres- 
sion fit assuming the relation 

ln[^(|TO,| + l)-£;(|TOf|)] =ifln(|m,| + l) + C, (32) 

where i^dm^l) is the energy with the one particle basis 
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1.031 




max(n) or max(R) 

FIG. 1: The figure shows the n-convergence of the two- 
particle ground state energy for a dot with a confining po- 
tential corresponding to fko = 3.32meV. The n-convergence 
is studied for a variety of starting points in the CCSD- 
calculation. The maximum n- value is varied from 1 to 
28, while |m;| < 1 for a Hartree-Fock starting point (red), 
different Local Density starting points (black), and a non- 
interacting starting point (blue with solid circles). The verti- 
cal lines with diamonds for max(n) — 9, 14 and 28 show the 
|m; |-convergence with max(|m;|) = [1, 10]. For reference we 
note that the energy using second order perturbation theory 
on top of Hartree Fock is 1.045 a.u.* when using the basis size 
(n, \mi\) < (28, 1). The (green) curve with hexagram markers 
shows the R = 2n+ |m; | -convergence when starting from the 
non-interacting basis set. 



TABLE IIL The |m;| -convergence (using all n, in this case 
n < 29) of the two electron ground state energy given in a.u.*. 
The extrapolated values are found according to the procedure 
described in the text. The extrapolated values agree well 
with the values 1.02164(1) and 1.02165(1) by Pederiva et al?? 
obtained through Variation Monte Carlo and Diffusion Monte 
Carlo methods respectively. For each extrapolated value three 
points were used in the linear fit. 



mi\ E{\mi\) \mi \ — )■ oo Upper bound Lower bound 



1 


1.024993 








2 


1.022767 








3 


1.022196 


1.021703 


1.022196 


1.021210 


4 


1.021971 


1.021667 


1.021971 


1.021363 


5 


1.021861 


1.021661 


1.021861 


1.021461 


6 


1.021799 


1.021660 


1.021800 


1.021520 


7 


1.021762 


1.021660 


1.021762 


1.021558 


8 


1.021737 


1.021661 


1.021738 


1.021584 


9 


1.021721 


1.021661 


1.021721 


1.021601 


10 


1.021709 


1.021662 


1.021709 


1.021614 



cut at maxdm^l) = \nie\ and K and C are the constants 
we find from the fit. The fits were made with three con- 



secutive differences at a time, which give the list of pre- 
dictions displayed in the third column in Table IIIII It 
is clear that this procedure can improve the results sub- 
stantially, especially when a rather low maximum \mi\ 
is used. The upper (lower) bound decreases (increases) 
monotonically as \mi\ increases and the extrapolated val- 
ues stabilize around 1.02166 a.u.*. This result agrees well 
with those obtained by Pederiva et a&^, 1.02164(1) and 
1.02165(1) a.u.*, through variational and diffusion quan- 
tum Monte Carlo methods respectively. 

Table IIVI shows the convergence of the ground state 
energies for two to eight confined electrons as functions 
of R = 2n -\- \mi\. This truncation scheme is a common 
choice in numerical studie o^^'^^ , and the motivation to 
study convergence as a function of this parameter is that 
the energy levels of a non-interacting two-electron dot 
are given by tnini = (2n + \mi\ + l)hbj. The confining 
potential corresponds to /iw = 3.32 meV, a strength that 
previously has been studied with Quantum Monte Carlo 
Methods^!, which are shown in the Table for compar- 
ison. A comparison with the two-electron results from 
Table IIIII indicates though that the radial convergence 
is substantially slower than the angular convergence, at 
least for this confinement strength. For more than two 
particles the angular convergence slows down due to the 
mutual repulsion felt by the electrons and the result- 
ing spread of the total wave function. This makes R 
a reasonable parameter for the truncation of the basis. 
We emphasize that we have used considerably higher R 
and/or more confined electrons than used in the available 
CI calculationsiSiil. The error relative to the Diffusion 
Monte Carlo method (expected to be the more accurate 
of the twc^^) is ~ 10~^ and on the same level or below 
as the difference between the two different Monte Carlo 
methods. For most purposes so far this would, we be- 
lieve, be a good enough convergence. We stress that the 
complete series of i? = [5, 15] for iV = 8 still did not 
take more than ^ 24 hours to compute on a standard 
desktop machine, which implies that further converged 
results can be obtained if necessary. 

However, it seems to be very difficult to improve the 
_R-cut results through extrapolation. All our attempts 
so far have resulted in quite unreliable predictions. For 
example, the extrapolated values did not stabilize, but 
showed a monotonously decreasing behavior. Instead we 
have investigated a third possible truncation scheme. In 
this scheme we have first used the Aitken's (5^-process, 
see e.g. Refi^, to speed up the convergence in n. After 
that we extrapolate to infinity in \mi\ as described at the 
beginning of this subsection in connection with Table Hill 
The results for the two-electron dot are displayed in Ta- 
ble|Vl where values for = 7, 8 and 9 for each |mi|— value 
were used to perform Aitken's convergence acceleration. 
In contrast to the attempted i?-extrapolation, the extrap- 
olated values do now stabilize, and the obtained value 
1.0217 a.u.* is very close to the result obtained in Ta- 
ble |llT] using a much larger basis set. The usage of this or 
similar extrapolation schemes for more than two confined 
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TABLE IV: 


The R = 2n + 




fnr the 9 — 


8 electron ground states for a confinin 


T^ni"pn1"ifil 

cL JJULCllLlCli 


pnvt'pQn(~i"n H i n o" fct 


h(jj — 3.32 meV. The results 


are given in a.u.* 












R 


2e" 


3e' 


4e~ 


oe 


6e 


7e" 


8e~ 


5 


1.02544 


2.24045 


3.72652 


o.oooi 1 


/ .DZD4U 


10.0913 


12.8065 


6 


1.02470 


2.23924 


3.72439 


o.o4oyu 


'T' 1 n /I o 


10.0624 


12.7284 


7 


1.02420 


2.23846 


3.72210 


0.04D/9 




10.0541 


12.7132 


8 


1.02383 


2.23791 


3.72207 


0.04404 


/ .OizU/ 


10.0498 


12.7067 


9 


1.02355 


2.23750 


3.72140 


0.04o/ 1 


/ .DiUUU 


10.0469 


12.7027 


10 


1.02333 


2.23719 


3.72089 


0.04234 


/ .dUo4d 


10.0448 


12.7000 


11 


1.02315 


2.23694 


3.72050 


^ ^41 fil 


1 .UU ( zo 


10.0431 


12.6979 


12 


1.02301 


2.23674 


3.72018 


5.54102 


7.60633 


10.0419 


12.6962 


13 


1.02289 


2.23657 


3.71991 


5.54054 


7.60557 


10.0408 


12.6949 


14 


1.02278 


2.23643 


3.71969 


5.54014 


7.60493 


10.0400 


12.6939 


15 


1.02269 


2.23631 


3.71950 


5.53979 


7.60438 


10.0393 


12.6930 


QMC 


1.02165(1) 


2.2395(1) 


3.7194(1) 


5.5448(1) 


7.6104(1) 


10.0499(1) 


12.7087(1) 


QMC'' 


1.02164(1) 


2.2339(1) 


3.7145(1) 


5.5338(1) 


7.6001(1) 


10.0342(1) 


12.6900(1) 



"Variation Monte Carlo, Pederiva et al^ 
'Diffusion Monte Carlo, Pederiva et al2^ 



TABLE V: The ground state energy of the hw = 3.32 meV 
two-electron dot obtained with an alternative extrapolation 
scheme: Aitken's (J'^-process is used on the n = 7, 8 and 9 val- 
ues for respective \mi \ to accelerate the n-convergence. Sub- 
sequently the |mi| — > cxD extrapolation is done as described in 
the text. The results are given in a.u.*. 

\mi\ E(|m;|) \mi\ ^ Qo Upper bound Lower bound 

1 1.025029" 

2 1.022809" 

3 1.022232" 

4 1.022001" 1.021697 1.021731 1.021658 

5 1.021893" 1.021703 1.021719 1.021687 

6 1.021839" 1.021734 1.021737 1.021697 

"Obtained by Aitken's 5^-process for n = 7 — 9. 

particles is an interesting topic for future studies. 



{hw K, 11.857nieV for GaAs material parameters) the 
results are for practical purposes exact when compar- 
ing with FCI and for A = 2 (^iw w 2.964meV for GaAs 
material parameters) the error is still never larger than 
^ 1.5 X 10^^ in units of hw. For > 6 the possibility 
to use much larger basis sets than in FCI-calculations 
is shown to be of uttermost importance. The errors in- 
troduced by truncating the basis sets in FCI-calculations 
are in many cases much larger than the error made by 
truncating some of the triple and quadruple excitations 
as done in CCSD. Moreover, when comparing with a Dif- 
fusion Monte Carlo study, for a potential strength close 
to what is estimated from the experiment by Tarucha et 
alji, the errors in the two to eight electron ground states 
are shown to be on the same level or less than the dif- 
ferences between Variational and Diffusion Monte Carlo 
results. This gives promise that the method, in future 
studies, can be used for the extraction of reliable infor- 
mation from experiments. 



V. CONCLUSIONS 

In conclusion, by comparison with results obtained 
with Full Configuration Interaction and two different 
Quantum Monte Carlo methods the Coupled Cluster Sin- 
gles and Doubles approach is shown to be a very power- 
ful method for two to eight electrons confined in a two- 
dimensional harmonic oscillator potential. For A < 1 
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